!
! Generates initial condition for running mhd2d code
! Needs mhd-ic.inp for parameters specification
! Also writes there some diagnostics
!
!=====================================================================

  MODULE size_params
    INTEGER ::   nx, ny
    INTEGER ::   nx2, nx1, nxx, ny2
    INTEGER ::   lx,ly,lxx
    INTEGER ::   n2d,n2dh 
    SAVE
  END MODULE size_params
!=====================================================================

  MODULE fftw_plans
!   For FFTw to work correctly, the size of these integers must
!   match the size of pointers in the language c.
      INTEGER, PARAMETER   ::   ikind = 8       !.. For alpha/64-bit machines
!     INTEGER, PARAMETER   ::   ikind = 4       !.. For most 32 bit machines
!     INTEGER, PARAMETER   ::   ikind = 7       !.. For g77 use (eg, PSR)
      INTEGER( kind=ikind) ::   plan1Df, plan1Db, plan2DRC, plan2DCR
!     INTEGER              ::   plan1Df, plan1Db, plan2DRC, plan2DCR
      SAVE
  END MODULE fftw_plans
!=====================================================================

     PROGRAM   ic

     USE size_params
     USE fftw_plans

     IMPLICIT NONE
     INCLUDE 'fftw_include.h'

     COMPLEX, ALLOCATABLE, DIMENSION (:,:)     ::  wk,wh
     REAL, ALLOCATABLE, DIMENSION (:,:)        ::  k2
     REAL, ALLOCATABLE, DIMENSION (:)          ::  kx,ky

     REAL :: time0, efinal, zk0v, qv, zk, zranv1, zranv2, wr, wi
     REAL :: ev, ens, pals, fact, gasdev
     INTEGER :: k_1, k_2, iseed, i, j, iwk


!------------------------------------------------------------

      OPEN(16,file='hd-ic.inp')
      read(16,*) nx
      read(16,*) ny
      read(16,*) k_1
      read(16,*) k_2
      read(16,*) efinal
      read(16,*) zk0v
      read(16,*) qv
      read(16,*) iseed
      read(16,*) iwk

      nx2=nx/2
      nx1=nx2+1
      nxx=nx+2
      ny2=ny/2

      lx=nx2
      ly=ny
!     lxx=lx
      lxx=lx+1

      ALLOCATE( wk(lxx,ny) )
      ALLOCATE( k2(lxx,ny), kx(lxx), ky(ny) )

      IF (iwk.eq.0) THEN

      ALLOCATE( wh(lxx,ny) )

        CALL set_plan2D (nx,ny)                ! FFTW plans
        n2d = ny * (nx + 2)
        n2dh = n2d / 2

      ENDIF

      CALL init_kays (kx, ky, k2)


      time0 = 0.0

      CALL srand(iseed)


!     Generate fourier coefficients wk

      wk = 0.0

      DO j=1, ny
      DO i=1, lxx
       
         zk = sqrt(k2(i,j))
         if ((zk.lt.float(k_1)).or.(zk.gt.float(k_2))) go to 1

        zranv1 = gasdev(iseed)
!... zran is gaussian with zero mean and unit variance
        wi = zranv1 * zk / sqrt(1. + (zk/zk0v)**qv )
        zranv2 = gasdev(iseed)
        wr = zranv2 * zk / sqrt(1. + (zk/zk0v)**qv )

        wk(i,j) = cmplx(wr,wi)

  1   CONTINUE      
      ENDDO
      ENDDO

!       Impose reality conditions
        CALL preal2D (wk)
        CALL zero_kyface (wk)
!       CALL zero_kxface (wk)


!... calculate the energy and globals
        CALL globals(wk,k2,ev,ens,pals)
        
!... normalize to efinal
        write(16,*)' '
        write(16,*)' COMPUTED INITIAL CONDITION CHARACTERISTICS '
        write(16,*)' unnormalized ev = ',ev
        fact  = sqrt(efinal/ev)


        wk = wk * fact

!...  calculate the energy and globals again
      CALL globals(wk,k2,ev,ens,pals)

      write(16,*)
      write(16,*)' AFTER NORMALIZATION'
      write(16,*) 'EKIN         =', ev
      write(16,*) 'ENSTROPHY    =', ens
      write(16,*) 'PALINSTROPHY =',pals
      close(16)
      

!     Writes initial condition
        IF (iwk.eq.0) THEN
        wh=wk
        CALL rfftwnd_f77_complex_to_real(plan2DCR,1,wh,1,n2dh,0,0,0)
        OPEN(1,file='w_0000',form='unformatted')
        CALL write_xarr (wh,time0)
        CLOSE(1)
        ELSE
        OPEN(1,file='wk_0000',form='unformatted')
        WRITE(1) time0,wk
        CLOSE(1)
        ENDIF


        END PROGRAM ic
!=====================================================================

      SUBROUTINE set_plan1D (n)

!     Set the FFTW 1D plans
!     Plans are passed to calling programs through the MODULE fftw_plans

      USE fftw_plans
      INCLUDE 'fftw_include.h'
      INTEGER, Intent( in)     ::  n

      CALL fftw_f77_create_plan( plan1Df, n, FFTW_forward,    &
                                FFTW_measure+FFTW_in_place)

      CALL fftw_f77_create_plan( plan1Db, n, FFTW_backward,   &
                                FFTW_measure+FFTW_in_place)


      END SUBROUTINE set_plan1D
!=====================================================================

      SUBROUTINE set_plan2D (n1,n2)

!     Set the FFTW 2D plans
!     Plans are passed to calling programs through the MODULE fftw_plans

      USE fftw_plans
      INCLUDE 'fftw_include.h'
      INTEGER, Intent( in)  ::  n1, n2

      CALL rfftw2d_f77_create_plan( plan2DRC, n1, n2,       &
                 FFTW_real_to_complex, FFTW_measure+FFTW_in_place)

      CALL rfftw2d_f77_create_plan( plan2DCR, n1, n2,       &
                 FFTW_complex_to_real, FFTW_measure+FFTW_in_place)


      END SUBROUTINE set_plan2D
!=====================================================================

      
      SUBROUTINE write_xarr (aa,time)

      USE size_params
      IMPLICIT NONE
      REAL, Intent(in)           ::   aa(2*nx1,ny)
      REAL, Intent(in)           ::   time
      INTEGER                    ::   i,j

      write(1) time,((aa(i,j),i=1,nx),j=1,ny)

      RETURN
      END SUBROUTINE write_xarr
!=====================================================================

      SUBROUTINE init_kays (kx, ky, k2)

!-------------------------------------------------------
!     kx = 0,1,2,...,nx/2  distributed on proc
!     ky = 0,1,2,...,ny/2,-(ny/2-1),...,-2,-1
!     k2 = kx*kx + ky*ky
!
!     Adjusted so that largest magn ky is POSITIVE
!-------------------------------------------------------

      USE size_params

      REAL, Intent( out), DIMENSION (lxx)        ::  kx
      REAL, Intent( out), DIMENSION (ny)         ::  ky
      REAL, Intent( out), DIMENSION (lxx,ny)     ::  k2

      INTEGER ::  i, j

      do i = 1, lxx
         kx (i) =  float(i-1)
      end do

      do j = 1, ny
         ky (j) =  float( MOD(j-1+ny/2,ny) - ny/2)
      end do
      ky (ny/2+1) =  float(ny/2)                       !.. Make ky_max > 0

      do j = 1, ny
        do i = 1, lxx
           k2 (i,j) = kx(i)*kx(i) + ky(j)*ky(j)
        end do
      end do

      END SUBROUTINE init_kays
!====================================================================

      SUBROUTINE psik (wk,fk,k2)

      USE size_params
      COMPLEX, Intent(in)      ::  wk(lxx,ny)
      COMPLEX, Intent(out)     ::  fk(lxx,ny)
      REAL, Intent(in)         ::  k2(lxx,ny)

      INTEGER                  ::  i, j

      do j=1,ny
         do i=1,lxx
            if (k2(i,j).ne.0.) fk(i,j)=wk(i,j)/k2(i,j)
         end do
      end do
      if (k2(1,1).eq.0.) fk(1,1)=0.0

      RETURN
      END SUBROUTINE psik
!====================================================================

      SUBROUTINE preal2D (aa)
!     for use with (kx,ky) MHD2D code
!
!     Imposes the interior AND edge reality conditions for a
!     pseudospectral algorithm.

      USE size_params
      IMPLICIT NONE

      COMPLEX,    Intent( inout)  :: aa(lxx,ny)
      INTEGER                     :: jy, jcjg
! ----------------------------------------------------

!     Modes (0, ky) are conj to (0, -ky)
      DO  jy = 2, ny/2
          jcjg = ny - jy + 2
          aa (1,jcjg) = CONJG ( aa (1,jy) )
      ENDDO

!     Do the special point k=(0,0) explicitly [=0, not just pure real]:
      aa(1,1) =  0.0

      RETURN
      END subroutine preal2D
!=====================================================================

      SUBROUTINE zero_kyface (aa)

      USE size_params
      IMPLICIT NONE

      COMPLEX,    Intent( inout)  :: aa(lxx,ny)

      aa (:, ny/2+1)  =  0.0

      RETURN
      END subroutine zero_kyface
!=====================================================================

      SUBROUTINE zero_kxface (aa)

      USE size_params
      IMPLICIT NONE

      COMPLEX,    Intent( inout)  :: aa(lxx,ny)

      aa (lxx, :)  =  0.0

      RETURN
      END subroutine zero_kxface
!=====================================================================

      SUBROUTINE globals(wk,k2,evt,enst,palst)

      USE size_params
      USE hd_params
      IMPLICIT NONE

      COMPLEX, Intent(in)      ::  wk(lxx,ny)
      REAL, Intent(in)         ::  k2(lxx,ny)
      REAL, Intent(out)        ::  evt,enst,palst

      INTEGER                  ::  i,j


      REAL :: wr, wi, ks

! evt   = kinetic energy = 1/2 sum (wk^2 / k^2)
! enst  = enstrophy = 1/2 sum (wk^2)
! palst = palinstrophy = 1/2 sum (k^2 * wk^2)


          evt   =  0.0
          enst  =  0.0
          palst =  0.0

          DO j=1,ny
             DO i=1,lxx

!        Count only the energy of the mode if ( kx > 0 .or. ky > 0 )
         if ( i.gt.1 .or. (j.gt.1 .and. j.le.ny2) ) then

                wr = real(wk(i,j))
                wi = aimag(wk(i,j))
                ks = k2(i,j)

                evt   =  evt   +  (wr*wr + wi*wi)/ks
                enst  =  enst  +  (wr*wr + wi*wi)
                palst =  palst +  (wr*wr + wi*wi)*ks

             ENDIF

         END DO
      END DO


      RETURN
      END SUBROUTINE globals
!=====================================================================

      FUNCTION GASDEV(IDUM)
      DATA ISET/0/
       common /gsave/gset
      IF (ISET.EQ.0) THEN
1       V1=2.*RAND()-1.
        V2=2.*RAND()-1.
        R=V1**2+V2**2
        IF(R.GE.1.)GO TO 1
        FAC=SQRT(-2.*LOG(R)/R)
        GSET=V1*FAC
        GASDEV=V2*FAC
        ISET=1
      ELSE
        GASDEV=GSET
        ISET=0
      ENDIF
      RETURN
      END
!=====================================================================



